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[57] ABSTRACT 

A computer-implemented method and apparatus for deter- 
mining position of a vehicle within 100 km autonomously 
from magnetic field measurements and attitude data without 
a priori knowledge of position. An inverted dipole solution 
of two possible position solutions for each measurement of 
magnetic field data are deterministically calculated by a 
program controlled processor solving the inverted first order 
spherical harmonic representation of the geomagnetic field 
for two unit position vectors 180 degrees apart and a vehicle 
distance from the center of the earth. Correction schemes 
such as a successive substitutions and a Newton-Raphson 
method are applied to each dipole. The two position solu- 
tions for each measurement are saved separately. Velocity 
vectors for the position solutions are calculated so that a total 
energy difference for each of the two resultant position paths 
is computed. The position path with the smaller absolute 
total energy difference is chosen as the true position path of 
the vehicle. 
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COMPUTER-IMPLEMENTED METHOD AND 

APPARATUS FOR AUTONOMOUS POSITION 
DETERMINATION USING MAGNETIC 
FIELD DATA 

5 

CROSS REFERENCE TO RELATED 
APPLICATIONS 

This application claims the benefit of priority under 35 
U.S.C. Section 119 (e) to U.S. Provisional application Ser. 

No. 60/035,949, filed Jan. 21, 1997. 

GOVERNMENT RIGHTS 

The invention described herein was made by an employee 
of the United States Government. It may be manufactured 15 
and used by or for the Government for governmental pur- 
poses without the payment of any royalty thereon or there- 
for. 

DOCUMENTS INCORPORATED BY 20 

REFERENCE 

This application is based in part on “Autonomous Navi- 
gation Recovery for Fine Pointing Low Earth Orbiting 
Spacecraft”, a dissertation submitted to the Faculty of the 
School of Engineering and Applied Science of The George 25 
Washington University in partial satisfaction of the require- 
ments for the degree of Doctor of Science, by Eleanor 
Ketchum Silverman, December, 1997, unpublished, and 
which is incorporated by reference. The following document 
is incorporated by reference, “Autonomous Deterministic 30 
Spacecraft Positions Using the Magnetic Field and Attitude 
Information” by Eleanor Ketchum, ION National Technical 
Meeting, Jan. 22-24, 1996, published in the publication 
“Proceedings of 1996 National Technical Meeting”, April, 
1996. The computer code in Appendix A is also incorporated 35 
by reference. “Autonomous Spacecraft Orbit Determination 
Using the Magnetic Field and Attitude Information” by 
Eleanor Ketchum, 19th Annual American Astronautical 
Society Guidance and Control Conference, Feb. 7-11, 1996, 
AAS 96-005 is also incorporated by reference. 40 

BACKGROUND 

The field of the invention is autonomous navigation 
systems for determining position based on inputs of mag- 45 
netic field and attitude data without the use of Kalman 
Filters. The present invention may be used for earthbound 
objects such as submarines and balloons as well as for 
spacecraft. 

In the context of autonomous spacecraft navigation, 50 
embodiments of the present invention are particularly useful 
in the context of low earth orbiting (LEO) spacecraft. All 
spacecraft navigation systems, at times, cannot determine 
their location. When this occurs, a spacecraft may have to 
rely on ground operations personnel to upload a position fix 55 
via a communications link. Because of the lack of navigation 
data, this communications link may be difficult for the 
spacecraft to establish. If the spacecraft cannot recover its 
orbit, instruments onboard the spacecraft may become 
damaged, or worse, the spacecraft maybe lost. To avoid this go 
situation, systems for a spacecraft to autonomously recover 
its position data, meaning without being supplied informa- 
tion externally, have been developed. 

One solution would be to put Global Positioning System 
Receivers on spacecraft. However, space qualified GPS 65 
Receivers have not as of the date of filing of this application 
been successfully demonstrated in space. In addition, the 


2 

power increase of fifteen Watts and the expense would be 
impractical for many spacecraft, particularly small LEO 
spacecraft. Moreover, many LEO spacecraft, such as those 
of the Small Explorer (SMEX) Program at the NASA 
Goddard Space Flight Center, only need to know their 
positions to within one hundred (100) kilometers, for con- 
tacting ground stations for instance, rather than to the 
accuracy in the hundreds of meters provided by GPS 
Receivers. 

Instead of putting new hardware on small spacecraft, 
software solutions used in conjunction with already existing 
hardware on spacecraft have been explored. One such 
method of determining position onboard a spacecraft used a 
transponder with either tracking satellite or ground station 
data. (See Gramling, C., et al, “TDRSS Onboard Navigation 
System (TONS) Flight Qualification Experiment,” Proceed- 
ings of the 1994 Flight Mechanics and Estimation Theory 
Symposium, Greenbelt, Md., May 1994, pp. 253-267 and 
Gramling, C. “Autonomous Navigation Integrated with 
NASA Communication Systems,” Proceedings of the 19th 
Annual AAS Guidance and Control Conference, 
Breckenridge, Colo., February, 1996, Paper No. AAS 
96-006.) Although the method is cost effective for transpon- 
der users, it can take 12 hours to a full day to converge on 
an initial position solution due to poor visibility. The method 
is also Kalman Filter based so that it requires a priori 
position information. 

Software methods to autonomously determine position, 
have been developed to use magnetic field data provided by 
magnetometers, which are already required by most space- 
craft for momentum management. Methods have been 
developed which sequentially filter magnetic field data pro- 
vided by magnetometers with extended Kalman Filters. 
Position estimates with the use of Kalman Filters have been 
shown to provide positions to better than fifty 50 km using 
flight data, (See Shorsi, G., and Bar-Itzhack, I., “Satellite 
Autonomous Navigation Based on Magnetic Field 
Measurements”, Journal of Guidance, Control, and 
Dynamics, Vol. 18, No, 4, July-August 1995, pp. 843-850) 
and even to better than one 1 km using simulated star tracker 
attitude information with simulated magnetic field measure- 
ments. (See Psiaki, M., “Autonomous Orbit and Magnetic 
Field Determination Using Magnetometer and Star Sensor 
Data,” August, 1993, AIAA-93-3825-CP.) 

However, these Kalman Filtering methods require an a 
priori knowledge of a position fix of the spacecraft for the 
filters to converge. Furthermore, the filters may take more 
than one hundred (100) minutes to converge on a solution 
and require extensive computing power. Computing power 
on a spacecraft must be shared with other necessary systems. 
For example, computing power onboard a SMEX spacecraft 
must be shared with every other subsystem, so it is safe to 
assume that only five 5 percent is available for the position 
determination process. In Psiaki, four iterations for a good 
guess and 25 iterations for a poor first guess (-2700 km 
error) which took 100 minutes, were obtained on a one 1 
Mflop workstation. An upper end SMEX only has about 900 
Kflops. 

Obtaining a position fix of within 100 km accuracy in 
significantly less time than 100 minutes based on magne- 
tometer data and attitude data without a priori knowledge of 
position is highly desirable in several contexts. For LEO 
spacecraft, the spacecraft would quickly be able to deter- 
mine its position and with that knowledge, recover its orbit. 
A quickly obtained position fix would provide a good first 
guess as the a priori knowledge of position needed by a 
Kalman Filter based method, thereby decreasing time to 
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convergence. Furthermore, a fine positioning spacecraft 
using GPS could also benefit from a real time coarse position 
estimate to decrease its time to a first position fix or offer a 
first when less than four (4) GPS space vehicles are visible. 

SUMMARY 

The present invention is a computer-implemented method 
and apparatus for finding a coarse position estimate for a 
vehicle using magnetic field measurements which addresses 
the above-mentioned shortcomings. The present invention is 
premised upon the Earth’s magnetic field being modeled in 
terms of nth order spherical harmonics. In particular, the 
predominant or core portion of the field can be represented 
by the gradient of a scalar potential function, which in turn 
can be represented by a series of nth order spherical har- 
monics. From this representation as the gradient of a scalar 
potential function, one may derive equations to solve for a 
magnetic field vector at a certain position in terms of that 
position in geocentric inertial coordinates. The relationship 
was inverted so that given a measurement of the magnetic 
field in inertial space and the time of its measurement, two 
position solutions can be deterministically found. To get a 
coarse estimate of the position, the dipole, or first order 
model of the geomagnetic field is used, hereinafter called the 
dipole characterization of the earth’s magnetic field. The two 
position solutions form an inverted dipole solution. 

An object of the present invention is to provide a 
computer-implemented apparatus for autonomous position 
determination of a vehicle comprising a memory having data 
and a program controlled processor having a clock for 
generating timing pulses. The data of the memory comprises 
measurements of magnetic field data for three axes mea- 
sured over a preset time period, times of the magnetic field 
measurements, and Gaussian coefficients of an nth order 
model of series of spherical harmonics representing the core 
portion of the magnetic field as a gradient of a scalar 
potential function. The processor comprises an inverted 
dipole solutions module which calculates an inverted dipole 
solution which comprises two position solutions. Each posi- 
tion solution is calculated for each measurement of the 
magnetic field and comprises a position vector for the 
vehicle and a vehicle distance from the center of the earth. 
Each inverted dipole solution is stored as data in memory. 
The processor further comprises a position determination 
module for determining a position for the vehicle from the 
data stored in memory. 

An object of the present invention is to provide the 
apparatus described further comprising a magnetometer for 
providing measurements of the magnetic field data and an 
analog to digital converter for converting the measurements 
of the magnetometer into computer readable form. The 
processor further comprises a translation module for obtain- 
ing magnetic field data in a fixed reference coordinate 
system such as earth centered, earth fixed or geocentric 
inertial coordinates by translating the computer readable 
magnetic data from body coordinates to such a fixed refer- 
ence coordinate system by using the attitude data. The 
translated magnetic field data is then stored in memory. 

Another object of the present invention is to provide a 
computer-implemented method for autonomous position 
determination of a vehicle comprising several steps. The first 
step comprises a program controlled processor inputting a 
measurement of three axes magnetic field data, said mag- 
netic field data being measured for a preset time period. 
Next, the processor determines whether there is a sufficient 
separation between the magnetic field vector of the mea- 


surement and a direction of a geomagnetic dipole according 
to a dipole characterization of the geomagnetic field in order 
for a correction scheme applied to a position solution to 
converge. The dipole characterization is a first order of a 
5 series of nth order spherical harmonics representing the 
geomagnetic field as a gradient of a scalar potential function. 

These two steps are repeated until it is determined that a 
sufficient separation exists. Upon a sufficient separation, said 
processor performs the next step of calculating an inverted 
10 dipole solution for each measurement of three axes magnetic 
field data. Each of the two position solutions of an inverted 
dipole solution are stored in a separate set in memory. These 
steps are repeated for each measurement of the magnetic 
field data measured for the preset time period thereby 
15 creating two sets of position solutions propagating two 
position paths for the vehicle with different directions. The 
processor then determines which of the position paths is the 
true position path of the vehicle. 

An object of the present invention is to provide a 
20 computer-implemented method for autonomous position 
determination further comprising these substeps, the first of 
which is determining two unit position vectors of the two 
position solutions in terms of the magnetic field vector and 
the direction of the dipole. The processor than calculates 
25 three possible solutions for the vehicle distance from the 
center of the earth from the magnetic field vector, the 
direction of the dipole and the total dipole strength of the 
geomagnetic field. Next the processor tests the three pos- 
sible solutions to obtain one having a real nonnegative value 
30 for the vehicle distance from the center of the earth. The 
processor then applies a correction scheme to each of the 
two position solutions of the inverted dipole solution until a 
stopping criteria is satisfied. 

35 Another object of the invention is to determine the true 
position path of the vehicle by the following steps. First, the 
processor determines a velocity vector for each position 
solution in each position path of the vehicle. Then a total 
energy for the position solution obtained from the first 
40 measurement of magnetic field data of the preset time period 
and a total energy for the position solution obtained from the 
last measurement of magnetic field data of the preset time 
period are calculated for each position path. A total energy 
difference is determined for each position path. The proces- 
45 sor then tests whether one total energy difference is greater 
than a noise factor, a constant predetermined value for the 
vehicle. Upon determining that one total energy difference is 
greater than the noise factor and one is less that the noise 
factor, the position path having the total energy difference 
50 less than the noise factor is selected as the true position path 
of the vehicle, 

Another object of the invention is to provide a successive 
substitutions method having a stopping criteria of conver- 
gence within a number of iterations as a correction scheme. 
55 Another object of the invention is to provide a computer- 
implemented apparatus for autonomous position determina- 
tion of a terrestial entity comprising a magnetometer for 
providing a measurement of magnetic field data for three 
axes, an analog to digital converter for converting measure- 
60 ments of magnetic field data into computer readable form, 
and a memory comprising a hemisphere location of the 
entity. The memory further comprising attitude data used by 
a processor comprising a translation module for translating 
the computer readable magnetic data from body coordinates 
65 to a fixed reference coordinate system. The translated mag- 
netic field data is stored in memory. Also stored in the 
memory are Gaussian coefficients of an nth order model of 
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series of spherical harmonics representing the core portion 
of the magnetic field as a gradient of a scalar potential 
function. The processor further comprises an inverted dipole 
solutions module for calculating an inverted dipole solution 
comprising two position solutions, each position solution 
comprising a position vector for the entity and an entity 
distance from the center of the earth, for each measurement 
of three axes magnetic field data. The processor also com- 
prises a position determination module for selecting a true 
position from the two position solutions based upon the 
hemisphere location stored in memory. 

Another object of the invention is to provide a computer- 
implemented method for autonomous position determina- 
tion of a terrestial entity comprising the following steps. The 
first step is obtaining a magnetic field measurement. The 
next is determining whether a sufficient separation exists 
between the magnetic field vector of the measurement and a 
direction of a geomagnetic dipole according to a dipole 
characterization of the geomagnetic field, to insure conver- 
gence of a correction scheme which will be applied to an 
inverted dipole position solution, wherein the dipole char- 
acterization is a first order of a series of nth order spherical 
harmonics representing the geomagnetic field as a gradient 
of a scalar potential function. Upon a determination of 
existence of the sufficient separation, the next step is to 
calculate an inverted dipole solution comprising two posi- 
tion solutions approximately 180 degrees apart, each posi- 
tion solution comprising a position vector for the entity and 
an entity distance from the center of the earth, for the 
measurement of three axes magnetic field data. Finally, a 
position is determined by selecting a true position from the 
two position solutions based upon a hemisphere location 
stored in a memory. 

Another object of the present invention is to provide a 
computer-implemented method for autonomous position 
determination of a terrestial entity comprising a magnetom- 
eter obtaining a magnetic field measurement in body coor- 
dinates. The measurement is then converted to computer 
readable form by an analog to digital converter. Next the 
magnetic field data is translated to a fixed reference coor- 
dinate system by using attitude data stored in a memory, said 
memory further comprising a hemisphere location of the 
entity. Next, it must be determined whether a sufficient 
separation exists between the magnetic field vector of the 
measurement and a direction of a geomagnetic dipole 
according to a dipole characterization of the geomagnetic 
field, to insure convergence of a correction scheme applied 
to a position solution, wherein the dipole characterization is 
a first order of a series of nth order spherical harmonics 
representing the geomagnetic field as a gradient of a scalar 
potential function. Upon a determination of non-existence of 
the sufficient separation, the previous steps are repeatedly 
performed until a sufficient separation exists between the 
magnetic field vector and direction of the dipole. Once the 
sufficient separation exists, the next step is calculating an 
inverted dipole solution comprising two position solutions 
approximately 180 degrees apart, each position solution 
comprising a position vector for the entity and an entity 
distance from the center of the earth, for the measurement of 
three axes magnetic field data. The last step is selecting a 
true position from the two position solutions based upon the 
hemisphere location stored in memory. 

Another object of the invention is to provide a successive 
substitutions method having a stopping criteria of conver- 
gence within a number of iterations being followed by 
applying a Newton-Raphson method for a second number of 
iterations as a correction scheme. 
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Further scope of applicability of the present invention will 
become apparent from the detailed description given here- 
inafter. However, it should be understood that the detailed 
description, while indicating preferred embodiments of the 
5 invention, is given by way of illustration only, since various 
changes and modifications within the spirit and scope of the 
invention will become apparent to those skilled in the art 
from this detailed description. Furthermore, all the math- 
ematical expressions are used as a short hand to express the 
10 inventive ideas clearly and are not limitative of the claimed 
invention. See Appendix B for a list of Nomenclature used 
in the mathematical expressions. 

BRIEF DESCRIPTION OF THE DRAWINGS 

15 These and other features, aspects, and advantages of the 
present invention will become better understood with regard 
to the following detailed description and accompanying 
drawings which are given by way of illustration only, and 
thus are not linmitative of the present invention, and 
20 wherein: 

FIG. 1A shows an embodiment of a computer- 
implemented apparatus for autonomous position determina- 
tion of a vehicle. 

25 FIG. IB shows an alternate embodiment of a computer- 
implemented apparatus for autonomous position determina- 
tion of a vehicle. 

FIG. 1C shows an embodiment of a computer- 
implemented apparatus for autonomous position determina- 
30 toion of a terrestrial entity. 

FIG. 2A shows a flowchart of a computer-implemented 
method for autonomous position determination of a vehicle. 

FIGS. 2B and 2C shows a flowchart of a computer- 
implemented method for calculating inverted dipole solu- 
35 tions. 

FIG. 2D shows a flowchart of a computer-implemented 
method for determining which of the position paths repre- 
sents the true path of the vehicle. 

40 FIG. 3 shows a preferred embodiment of a computer- 
implemented method for autonomous position determina- 
tion of a vehicle. 

FIG. 4 shows an embodiment of a computer-implemented 
method for autonomous position determination of a terres- 
45 trial entity. 

DETAILED DESCRIPTION 

FIG. 1A shows an embodiment of an apparatus for 
autonomous position determination of a vehicle. The appa- 
50 ratus (12) consists of a program controlled processor (22) 
and a memory (14) storing the data needed by and generated 
by the software program embedded in the processor for 
determining the position of a vehicle from magnetic field 
data. The processor has an Inverted Dipole Solutions Mod- 
55 ule (24) and a Position Determination Module (32), which 
are embodied in software. The memory (14) stores three 
axes magnetic field data (16), attitude data (30), measure- 
ment times of magnetic field data (18), and Gaussian coef- 
ficients for an nth order spherical harmonics model of the 
60 geomagnetic field. 

FIG. IB is a preferred embodiment of the invention in 
which the magnetic field data (16) is measured by a three 
axes magnetometer (34). An analog to digital converter (36) 
converts the magnetometer measurements to a computer 
65 readable form. A translation module (38) in the program 
controlled processor translates the magnetometer data from 
body coordinates to geocentric inertial coordinates using the 
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attitude data. FIG. 2 shows a computer implemented method 
for obtaining inverted dipole solutions, each of which com- 
prises two position solutions 180 degrees apart. The method 
steps may be easily coded into a software program. Before 
describing the steps of the method, a discussion of their 5 
development follows to understand how the method works. 

As stated in the Summary, the magnetic field of the Earth 
can be represented by the gradient of a scalar potential 
function. The function can therefore be modeled by a series 
of k th order spherical harmonics: 1C 
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Where a is the equatorial radius of the earth, r, 0, and (|) are 
the geocentric, co-elevation, and longitude from Greenwich 
for any point in space. From this model, the magnetic field 
vector at a position may be defined in terms of the position. 20 
An nth order field model that may be used is the 10th order 
model incorporating the 1990 International Geomagnetic 
Reference Field Gaussian coefficients, g n m and h n m . The 
Lengendre functions, P m n , are assumed Schmidt normalized 
(see Plett, M., “Magnetic Field Models” Spacecraft Attitude 25 
and Control , J. R. Wertz (editor), Kluwer Academic 
Publishers, Dordrecht, The Netherlands, 1978, pp. 779-786, 
and Chapman, S, and Bartels, J., Geomagnetism, Oxford 
England, Clarendon Press, 1940) such that the P M m satisfy: 


Therefore, given a position in geocentric coordinates, it is 
possible to compute the magnetic field vector at that position 
(see Plett): 
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( 2 ) 

or in geocentric inertial coordinates: 


where 8 is the Kronecker delta. Further, the Gaussian 35 
normalized Legendre functions, pn ’ m , are related to the 
Schmidt normalized Legendre functions by: 

P n m =S n , m ? n ’ m 

where 

„ [(l-<5»)(«-m)!l 1/2 (2«-l)!! 

(B + m)! J (n-m)< 

where (2n-l)!! =1*3*5 . . . (2n-l). Plett points out that these 
factors, S nm , are best combined with the Gaussian coeffi- 
cients since they are independent of r, <|), and 0 and therefore 
need only be computed once. Therefore the following rela- 
tionships can be defined: 


(3) 

40 

(4) 

45 


B^B,. cos 6+B 0 sin 6)cos a-B^ sin a 

(18) 

B^B,. cos 6+B 0 sin 6)sin a+B^ cos a 

(19) 

B z =(B f sin 6-B 0 cos 6)cos a 

(20) 


Where 8 is the declination and a is the right ascension of the 
vehicle, where the right ascension is the sum of the longitude 
and the right ascension of the Greenwich meridian. 

These equations may be inverted so that, given a mea- 
surement of the magnetic field in inertial space, a position 
can be found. To get a coarse estimate of the position, the 
dipole model for the earth’s magnetic field was considered. 
The dipole characterization of the earth’s magnetic field in 
geocentric inertial coordinates is: 


B x = 


a 3 H 0 
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[3(m ■ R)R X - sin0 OT cosar OT ] 


(21) 


iT’ m =s n , m g n m (5) 

h n ’ m =S n ,Jln m ( 6 ) 

In order to compute these functions, the following recursive 
relationships are useful (see Plett): 
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(25) 

and 
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(26) 
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-continued 


a 3 Hq = a 3 [, g J 2 + g\ 2 +h\ 2 ] / = the total dipole strength 


(27) 


Where R ^ are the unit directions of the vehicle position, 
r is the vehicle distance from the center of the earth, R is the 
unit direction of the vehicle, a G0 is the right ascension of the 
Greenwich meridian at the time from which t is measured, 

da G 

~df 


is the average rotation rate of the earth (360.9856469 
deg/day), and 0 m are the east longitude and co-elevation 15 
of the dipole, and a m is the right ascension of the dipolw. 
The therefore m represents the direction cosine of the dipole. 

Given B and the location of the dipole at the desired time, 
the quadratic equations (21-27) can be solved for the 
direction cosines: 20 

2\B\m ™ 

— — — +B 

R = ± ' ' T 25 

{ m { m A ->)}* 

= — +m-B 

{(m-B) + ^D {(m-B)+^D ) / 


and the vehicle distance from the center of the earth, r: 30 


r 3 = 


a 3 H 0 

s(B-B) 


(m-B + ^Jd") 


(29) 


where: 


35 


D=(m-B) 2 +8|B| 2 (30) 

Even though equation (29) is 3rd order, only a single real 40 
non-negative solution exists. 

These expressions (28-30) form the “inverted dipole” 
solutions. The ± in the expression for the direction cosine 
solution in expression 28 implies that there are 2 position 
vector solutions to these equations 180° apart. Further 45 
information is used to determine the correct solution. 
Furthermore, a singularity exists mathematically when the 
magnetic field vector is 180° from m, the direction of the 
magnetic dipole. 

These inverted dipole solutions can be in error by more 50 
than 1000 km if the true position is near the magnetic 
equator and 100 km elsewhere. This error is due to the fact 
that the higher order components have a large effect near the 
equator. The dipole model is first order so that the higher 
order terms have been neglected thus far. As a consequence 55 
of the contribution of the higher order terms, particularly 
near the magnetic equator, a correction scheme to the 
inverted dipole solution was developed. 

FIG. 2A shows a flowchart of the steps to be performed 
by an embodiment of the computer-implemented method for 60 
autonomous position determination (200). A processor 
inputs(208) measurements (206) of the three axes magnetic 
field data in geocentric inertial coordinates. These measure- 
ments are no older than the beginning of a preset time 
period. The preset time period is discussed below in the 65 
discussion on velocity vectors. Next a determination (210) 
as to whether a sufficient separation exists between the 


magnetic field vector and the dipole is made. If a sufficient 
separation exists, the processor calculates an inverted dipole 
solution for the measurement (212). 

FIGS. 2B and 2C outline steps that a program controlled 
processor performs to calculate the inverted dipole solution. 
The direction of the dipole m is calculated (222). Next two 
unit position vectors are determined from the magnetic field 
vector of the measurement and m (224). Equations 28 and 30 
may be embedded in software instructions to determine 
these unit position vectors, which will be used by the vehicle 
to determine its position. Next, three possible solutions for 
the vehicle distance from the center of the earth are calcu- 
lated by the processor (226). Equation 29 may be pro- 
grammed in software to calculate three possible solutions for 
the vehicle distance from the center of the earth. The total 
dipole strength may be treated as a constant term with value 
7.84e+015 Wb m, although the dipole strength has 
decreased by about 1.28% from 1975 to 1990. Next, the 
three solutions are tested to see which represents a real 
physical distance by determining which of the three solu- 
tions has a real nonnegative value (228). Each unit position 
vector and the vehicle distance represent a possible position 
solution, so that the inverted dipole has two position solu- 
tions. For the reasons given above, a correction scheme is 
applied to each position solution until a stopping criteria is 
satisfied (230). 

Both a successive substitutions method and a successive 
substitutions method in combination with a Newton- 
Raphson method are used in embodiments of the invention. 
The preferred embodiment of FIG. 3 processes each inverted 
dipole solution through a stopping criteria of twenty itera- 
tions for the successive substitutions method and five itera- 
tions for the Newton-Raphson method. The Newton- 
Raphson method is known for its quadratic convergence 
properties but is not robust to poor initial estimates. 
Alternatively, successive substitutions is known for robust- 
ness to errors in initial estimates, but has linear convergence 
properties (see Battin, R. A., An Introduction to the Math- 
ematics of Methods of Astrodynamics, New York, American 
Institute of Aeronautics and Astronautics, Inc., 1987, pp. 
192-217). 

Successive substitutions is an iterative solution. The mag- 
netic field vector is divided into the dipole component and 
higher order terms (HOTs), or: 

B > 1 =DIPOLE+HOT (31) 

and consequently, 


HOT= B j-DIPOLE 
Where: 


(32) 


B j-A B BCS (33) 

and B bcs is the magnetic field vector measured in body 
coordinates and A is the direction cosine matrix of the 
attitude data which transforms from body coordinates to 

inertial space. The quantity B a is measured. Onboard a 
spacecraft, for instance, the attitude direction cosine matrix, 
A, is determined using measurements from a star tracker, 
sun sensor or other sensor whose measurements may be used 
to determine attitude. The present invention only requires 
attitude data as an input and is not dependent on the source 
from which is was obtained. The magnetic field vector is 
measured by a three axis magnetometer. 
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The dipole inversion logic above yields the new operator, 
DIPOLE _1 ( B ), which can be used iteratively to improve R : 

R^^DIPOLE-^Bi-HOTCR^^R,) (34) 

5 

To get the HOT(R^), the R* estimate is used to compute 
both the nth order field and the dipole field in equation (32). 

In a preferred embodiment of the invention, the nth order 
model is the 1990 10th order model. 

In order to examine the convergence behavior of this 10 

iterative scheme, R is rotated into magnetic coordinates so 
that m is now e 3 and: 


- a (a 3 H 0 )i 

] 


^ T 

B\B\ L/5+8|5| +B 


12\B\ +6B . 


B + 8|B| + B 


i = x, y 


- a (a 3 H 0 ) 

^)=(— j 


3 A ” 3 - 

\B\ U/f + 8|5| +B 


1^> 2 l^Z 

-\B\ +-B 

3 6 


B + 8|5| + B 

2 


(35) 


(36) 


15 


20 


25 


where B is the measured magnetic field vector in inertial 30 
space minus the computed higher order terms. 

Clearly, equation (36) shows that a singularity exists when 
l\=- L. The following analysis will provide insight into how 
close to -1 B, can be and still have guaranteed convergence 
for the successive substitutions scheme. A sufficient condi- 35 
tion for successive substitutions of a system of n non-linear 
equations to converge is that the absolute value of each of 
the partial derivatives of F(R) in equation (34) with respect 
to each variable is less than 1/n, where n is the number of 
variables. (Carnahan, B., Luther H. A., and Wilkes, J. O., 40 
Applied Numerical Methods, New York, 1969.) In the case 
where n=3, results from these calculations show that the 
sufficient condition for successive substitutions convergence 
is met if the measured unit magnetic field z direction in 
magnetic coordinates is greater than -0.978, or 12° from the 45 
magnetic dipole direction. Near the magnetic equator, the 
magnetic field direction changes at a rate of approximately 
2.907° latitude. Therefore, the dipole changes 24° with a 


where J is the 3x3 Jacobian, 
dB 

dl' 


which can be derived analytically (see Appendix C). 

The Jacobian has properties which are useful from a 
computation aspect. The matrix is symmetric and the matrix 
has zero trace. Therefore, only five elements need to be 
computed to construct the Jacobian. 

While the successive substitutions requires the derivative 
of the function to be less than 1/n to guarantee convergence 
for a system of n equations, the Newton-Raphson method 
simply requires the measured magnetic field vector be 
sufficiently separated from the dipole direction to allow for 
measurement uncertainty. Therefore, if the sufficient condi- 
tion for the successive substitutions is met, the Newton- 
Raphson method’s requirement that the measured magnetic 
field vector be sufficiently separated from the dipole is met. 

Near the magnetic equator, the field is at a minimum. In 
low earth orbit, the field minimum is approximately 225 
mgauss. If measurement uncertainty is 0.7 mgauss per axis 
in the magnetometer, or 1.2 mgauss RSS, the field vector 
must be greater than atan(l .2/225) or 0.3° from the negative 
dipole direction. Further, there is approximately a one 
degree error in the spherical harmonic model at any given 
point, due to non core magnetic activities. Therefore, the 
measured magnetic field direction must be at least 1.3° from 
the negative dipole direction. Evaluating the Jacobian in 
Appendix C at the magnetic dipole equator for a 1° change 
in the spacecraft latitude yields a maximum 3.9° change in 
the magnetic field. 

The relationship between the error in the measurement of 

the magnetic field, A B B csmeasured> as well as the error in the 
spacecraft attitude, AA, to the resulting error in the position 
solution is developed as follows: 


0=F( R+AR)=B 1 modeled" (A+AA) B BCSmeasured- (A+ A) A 

V BCSmeasured -KOT (39) 

or, neglecting the higher order terms: 

B 1 modeled^ ^ R ~(A+AA) B BCSm easured ± (A+A A)A B BCSmeasuredAfy 

which reduces to: 

A R «J 1 (AA B BCSmeasurec&A* * ^ B BCSmeasured ) (41) 


change of approximately 8° in latitude. In sum, R x and R y 5Q 
both have singularities as B^ approached -1. However, as 
long as B^ is greater than -0.979, or 12° from the magnetic 
dipole, successive substitutions will converge on a position 
solution. 

The Newton-Raphson corrector, while more sensitive to 55 
poor initial conditions than successive substitutions, gener- 
ally converges more rapidly (see Battin). The analytical 
partial derivatives required for the method can be computed 
in a straight forward way. 

For the Newton method, the function that must go to zero 6Q 
is defined as: 

F(R)=b‘ 1 modeled(R)- B x (37) 


Therefore, the errors in R are directly proportional to J -1 . 
When the gradient in the field is small, the error is large. So, 
the position error is less effected by measurement errors at 
higher latitudes. Assuming the error in attitude is not sig- 
nificant (i.e. is less than 100 arcsecs RSS for a single star 
tracker scenario), 

AR=r 1 (AAB) (42) 

Let \ min be the minimum absolute valued eigenvalue of J. 
Then a singularity occurs in the limit as \ min goes to zero, 
where J becomes non-singular. A more easily computed 
indicator that J is becoming non-singular is the determinant 
of J approaching zero. The matrix J can be decomposed as 
follows: 


Again, B 1 is measured. Then the Newton iterative scheme 
is: 


R* +1 =R*-r 1 F(R A ) 


J=M -1 DM (43) 

where D is a diagonal matrix in which the eigenvalues of J 
are on the diagonal, and the columns of M are the eigen- 


(38) 
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vectors of J. Substituting equation (43) into equation (42) 
and multiplying both sides by M yields: 

MAR=D -1 M(AAB) (44) 

Multiplying together the elements of the vectors on each 
side of equation (44) yields: 

(AR 1 -AR 2 -AR 3 )=(AB 1 -AB 2 -AB 3 )/|J| (45) 


where R, and B, are simply transformed into the eigenspace. 

For example, to assure that the position solutions have errors 
less than say 125 km in each direction, and if the measure- 10 
ment uncertainty is 6 mgauss in each direction, then the 
determinant of J must be >0.0001. Further, equation (42) 
represents a worst case bound on the error since the higher 
order terms were neglected. 

After the inverted dipole solution has been calculated and 15 
corrected (212), each position solution is stored in a separate 
data set in memory (214) so that, once it is determined that 
there are no more measurements for the preset time period 
(216), two position paths are stored in memory. Of these two 
position paths, the true position path of the vehicle must be 
determined in order for the vehicle to know where it is (218) 
so that it can act appropriately. 

For terrestial entities including living and nonliving 
things, with knowledge of which hemisphere the entity is 
located, one of the position solutions from the inverted ^ 
dipole solution can be selected as the true position on one 
measurement. FIG. 1C shows a computer-implemented 
apparatus for autonomous position determination of a ter- 
restrial entity. The memory (14) comprises the hemisphere 
location (35) as well as Attitude Data (30), the Gaussian ^ 
coefficients(20) and a translated magnetic field measurement 
(25). FIG. 4 shows the flowchart of a computer-implemented 
method for autonomously determining the position of a 
terrestrial entity. Step (402) is obtaining a magnetic field 
measurement. Step (404) is determining whether a sufficient 35 
separation exists between the magnetic field vector of the 
measurement and a direction of a geomagnetic dipole 
according to a dipole characterization of the geomagnetic 
field, to insure convergence of a correction scheme to an 
inverted dipole position solution, wherein the dipole char- 
acterization is a first order of a series of nth order spherical 
harmonics representing the geomagnetic field as a gradient 
of a scalar potential function. Step (406) calculates an 
inverted dipole solution which combined with the hemi- 
sphere location in Step (408) allows selection of the true ^ 
position. 

In order to determine the true path, one method takes 
advantage of the fact that the magnetic field is not sym- 
metrical about its equator so that the incorrect or ambiguous 
position path will not follow the same Keplerian orbit that 
the true path does. Further, it is assumed to be unlikely that 
the ambiguous solution will follow a Keplerian orbit at all. 
The total energy of an orbit can be computed using just the 
magnitude of the position and velocity vectors as follows: 


In order to calculate a total energy, velocity vectors are 
required. In FIG. 2D, a velocity vector for each position 
solution in each of the two position paths is determined 
(232). For a spacecraft, one way to calculate velocity vectors 
is by using the Taylor series expansion of the equations of 
motion, the velocity can be computed from at least two 
position vectors. 

Battin describes a method for propagating an orbit from 
an initial position and velocity vector. Since the equations of 
motion for the two body problem generally do not have a 
closed form analytic solution, the Taylor series expansion is 
considered. The position of a body can be expressed as: 

R(f)=F(f)R 0 +G(f)V 0 (47) 

where the coefficients F(t) and G(t) can be represented as a 
power series in (t-t 0 ). The coefficients each satisfy the 
differential equations: 


d 2 F(t) 
dt 2 

d 2 G(t) 
~dfi~ 


+ sF{t) = 0 
+ sG(t) = 0 


(48) 

(49) 


Along with the initial conditions that F(t0)=l, G(t0)=0, 


dF 
dt tQ 


dG 

- 0, and 

dt 


lf o 


= 0, 


yields: 


F(t) = Y J Fn(i-rof 

n= 0 


(50) 


G(t) = Y j G„(t-t 0 r 

n= 0 


(51) 


These coefficients are given by: 

F 0 = 1 Fi = 0 F 2 = -_£ 0 F 3 = _£ 0 A 0 
F 4 = -5 / 8 £ 0 Ao + 1 / 8 £ 0 |Ao - 1 / 12«o etc. 

Go = 0 G\ = 1 G 2 = 0 G 3 = — 1 / 6 f>o G 4 = — 1 / 4f>oAo etc. 


where e, X, and ip are the Legrange fundamental invariants: 


£ = 


r* 


(52) 



(53) 


TE = 



55 


(46) 



(54) 


where v is the magnitude of the velocity, r is the magnitude 
of the position vector, and fd is the earth’s gravitational 
constant. The difference of the total energy at time t and time 
t+n minutes, ATE, should be equal to zero plus some noise 
due to the measurement errors. Therefore, to break the 
ambiguity the solution with the smaller abs(ATE) should be 
chosen. However, further quantification is required for the 
cases where either both abs(ATE) solutions are large or both 
close to zero. 


For this application, spacecraft position vectors can be 
6Q directly computed using the method described above, but 
there is no independent measure of velocity. An estimate of 
the spacecraft velocity can be computed in a two step 

process. First, equation (47) can be solved for V 0 given at 
least two position fixes and an initial estimate of V 0 to 
65 compute the approximate F and G series. Further, with more 
than 2 position fixes, a least squares fit can provide improved 
velocity vectors: 



15 
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V 0 =(R(t)-F(t)R 0 )/G(t) (55) 

In the particular embodiment of FIG. 3 using the equations 
of (47) through (55) in embedded software in the processor, 5 
the F and G series are intentionally at 3 rd order to eliminate 

Rand consequently minimize the occurrence of V on the 
right hand side of equation (55). Also note that in equation 
(55), the coefficients F(t) and G(t) are actually nxl matrices 10 

where n is the number of R’s used to solve for V 0 . The 
vector divide in equation (55) is more precisely defined as 
multiplying by the pseudo inverse of G(t). This becomes a 
least squares solution to the system of 3xn equations with 3 
unknowns. The second step is to include more powers of 15 

(t-t 0 ) by iteratively using the computed V 0 to determine 
new Legrange invariants. Subsequently, equation (55) can 
be reevaluated. 

The maximum amount of time over which R vectors 2Q 

should be saved and used to compute V 0 can be quantified 
as follows. The truncated F and G series is a poor estimator 
of the orbit if t-t 0 becomes too large. Since the series are 
truncated Taylor series, the nth term is proportional to 
(t-t 0 ) n . Since the A th term is being truncated here, it is 25 
necessary to examine where that term’s influence will be 
significant with respect to the desired error tolerance. The 
present invention is attempting to find positions on the order 
of 100 km error. A significant influence of the (t-t 0 ) 4 term 
would be considerably less than 100 km. A reasonable limit 30 
is 20 km, so that it is necessary to define (t-t 0 ) such that: 


F 4 (t-t 0 ) 4 R 0 +G 4 (r-r 0 ) 4 V o >20 km (56) 

In low earth orbit, F 4 is on the order of 10 -14 , R 0 is 7000 km, 

G 4 is on the order of 10 -13 , and V 0 «7 km/sec leads to (t-t 0 ) 35 
on the order 700 seconds. Therefore, the maximum amount 
of time over which position vectors should be saved and 
used to compute velocity vectors using this method is 
approximately 10 minutes. 

The minimum time between two position vectors used to 40 
compute velocity is related to the velocity error requirement. 
When using only two position vectors, and without iteration 
equation (55) reduces to the difference in the vectors to 

determine V 0 . Therefore, the worst case error in velocity is 45 
the error in the difference in the two positions divided by the 
time between them. These values are limited by the desired 
performance. Further, when more than two position vectors 

are used in equation (55), V 0 is determined by a least 
squares approximation and the error is reduced as a function 50 
of 1/^n where n is the number of position vectors used. This 

is a conservative estimate as the errors in B are correlated. 

The embodiment of the invention in FIG. 3 (300) collects 
five minutes magnetometer data (302) from which position __ 
vectors were calculated and will be used to determine 
velocity. This value is less than the 10 minute recommended 
maximum value for the preset time period. 

With velocity vectors determined, a magnitude for veloc- 
ity can be determined for computing the total energy of each 
position path. 



Again, v is the magnitude of the velocity, r is the magnitude 
of the position vector, and // is the earth’s gravitational 
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constant. The difference of the total energy at time t and time 
t+n minutes, ATE, should be equal to zero plus some noise 
due to the measurement errors. Therefore, to select the true 
position path, the position path with the smaller abs(ATE) 
should be chosen. However, further quantification is 
required for the cases where either both abs(ATE) solutions 
are large or both close to zero. Assuming that there are errors 
in measurement, the difference between the total energy at 
two different times, t 0 and t„, in orbit is: 




V'O +r '0 Sr 'n +r 'nK 

( 8 r tn - Sr t ) 

(<5v, 0 v fQ - 8v t v t ) + fi\ — 

\V’0 +r '0 S % + r ’nK) 


■■ TE 'o ~ TE 'n + 2 (|5V '0 V '0 “ Sv <n V <n) H 
Sr < 0~ Sr 'n 


V<0 +r <0 Sr <n + r ’n Sr >o 


(57) 


Therefore if: 

ATE=TE, 0 -TE, n +6TE 

Then combining equations (57) and (58) yields: 


( 8r t - 6r t 

STE = {8v t v, - Sv t v t ) + fi — 

K'Vo +^A +r Aj 


(58) 


(59) 


Since it has been assumed that the spacecraft is in low 
earth orbit, any orbits must be nearly circular. The maximum 
eccentricity possible below 1000 km altitude has been 
shown to be 0.06. For a circular orbit, v and r are constant 
and 6v and 8r are latitude dependent. Worst case values for 
a given orbit for a given 8v and 5r can be estimated and 5TE 
calculated once. The quantity 8r ?o -8r^ should equal the 
maximum 5r*2 since that is the worst case difference 
between the two assuming the error in r is maximum positive 
at the end of the span and maximum negative at the 
beginning. Similarly, should equal 8v*v*2. 

Therefore, in the worst case, equation (59) becomes: 


8TE = 2<5vv + 



(60) 


The true positon path cannot be definitively resolved until 
one of the abs(ATE) is greater than the error caused by the 
measurement noise as defined in equation (60). 

Defn. A definitive ambiguity resolution occurs when 
exactly one solution is less than 5TE. Further, if both 
abs(ATE) values are greater than 5TE, then the ambi- 
guity cannot be definitively broken. 

For example, for a spacecraft in a circular orbit with an 
inclination of 82°, the maximum change in latitude a space- 
craft needs to traverse for a definitive ambiguity resolution 
is 40°. However, the mean change in latitude is 4°. For that 
particular orbit, that change in latitude translates to a mean 
time of one minute. Simulations will also show that the 
instance of an incorrect choice of solution when one abs 
(ATE) is less than 5TE and the other is greater than 5TE 
occurs only within 10° of the magnetic equator. 

A constant value can be used for 5TE based on known 
error measurements for a vehicle so that 5TE is a predeter- 
mined noise factor. 

In step (234) FIG. 2D, the program controlled processor 
calculates a total energy for the position solution obtained 
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from the first measurement of the magnetic field data of the 
preset time period, and in step (236) another total energy for 
the last position solution of the preset time period. Step 
(238) computes a total energy difference for each position 
path. When both total energy differences have been com- 5 
puted (240), the differences are compared with a predeter- 
mined noise factor (242). If one is greater than the noise 
factor, that is not the true position path, so the other one is 
selected as the true position path (244). If neither difference 
is greater than the noise factor, than a true position path 
cannot be selected. The program controlled processor inputs 
another measurement and repeats the process (208). 

Each position and velocity vector can be propagated to 
predict the positions up to about 10 minutes into the future 
using the above scheme, or even further as more (t-t 0 ) terms 
are included. This scheme includes only the forces of the two 15 
bodies. The largest additional disturbing force is J 2 , or the 
force due to the oblateness of the earth. Over an orbit 
average, J 2 tends to precess the ascending node as well as the 
argument of perigee (see Wertz, J. R. “Summary of Orbit 
Properties and Terminology,” Spacecraft Attitude Determi- 20 
nation and Control , J. R. Wertz (editor), Kluwer Academic 
Publishers, Dordrecht, The Netherlands, 1978, pp. 67-69). 

In order to study the short term effects of J 2 , a Monte Carlo 
approach was taken. For 1000 randomly chosen initial 
conditions, low earth orbits were propagated ten minutes 2 s 
ahead using two different propagators. One propagator 
included only two body effects, while the other included J 2 . 
The differences at the end of the ten minute spans were then 
examined statistically. The mean RSS difference was 
approximately 7 km, the standard deviation was 2 km, and 
the maximum was 13 km. When considering errors on the 
order of less than 50 km, a simple propagator with J 2 should 
be used. Therefore, during periods of time when the space- 
craft passes near the magnetic equator and no deterministic 
solution is available, the orbit can be propagated. 35 

Obviously, many modifications and variations of the 
present invention are possible in light of the above teach- 
ings. For example, the methods are not dependent on a 
particular coordinate system so different coordinate systems 
may be used such as geocentric inertial coordinates or earth 
centered, earth fixed coordinates. It is therefore to be under- 
stood that within the scope of the appended claims, the 
invention may be practiced otherwise than as specifically 
described. 

What is claimed is: 

45 

1. A computer-implemented apparatus for autonomous 
position determination of a vehicle comprising: 

a memory having data, said data comprising measure- 
ments of magnetic field data for three axes represented 
in a fixed reference coordinate system, times of the 5Q 
magnetic field measurements, and Gaussian coeffi- 
cients of an nth order model of series of spherical 
harmonics representing the core portion of the mag- 
netic field as a gradient of a scalar potential function; 
said memory being accessible by a program controlled 55 
processor having a clock for generating timing pulses; 
said processor comprising an inverted dipole solutions 
module for calculating an inverted dipole solution 
comprising two position solutions, each position solu- 
tion comprising a position vector for the vehicle and a 60 
vehicle distance from the center of the earth, for each 
measurement of three axes magnetic field data; 
said inverted dipole solution being stored as data in said 
memory; and 

said processor further comprising a position determina- 65 
tion module for determining a position from the data 
stored in said memory. 
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2. The apparatus of claim 1 wherein said fixed reference 
coordinate system is geocentric inertial coordinates. 

3. The apparatus of claim 1 wherein said fixed reference 
coordinate system is earth centered, earth fixed coordinates. 

4. The apparatus of claim 1 wherein the Gaussian coef- 
ficients comprise all the coefficients of the 1990 10th order 
model of the International Geomagnetic Reference Field. 

5. The apparatus of claim 1 further comprising: 

a magnetometer for providing measurements of magnetic 
field data for three axes; 

an analog to digital converter for converting measure- 
ments of magnetic field data into computer readable 
form; 

said memory further comprising attitude data; 

said processor further comprising a translation module 
translating the computer readable magnetic data from 
body coordinates to a fixed reference coordinate system 
by using the attitude data; and 

said translated magnetic field data being stored as data in 
memory. 

6 . The apparatus of claim 5 wherein the Gaussian coef- 
ficients comprise all the coefficients of the 1990 10th order 
model of the International Geomagnetic Reference Field. 

7. The apparatus of claim 5 wherein said fixed reference 
coordinate system is geocentric inertial coordinates. 

8 . The apparatus of claim 5 wherein said fixed reference 
coordinate system is earth centered, earth fixed coordinates. 

9. The apparatus of claim 7 wherein the Gaussian coef- 
ficients comprise all the coefficients of the 1990 10th order 
model of the International Geomagnetic Reference Field. 

10. The apparatus of claim 8 wherein the Gaussian 
coefficients comprise all the coefficients of the 1990 10th 
order model of the International Geomagnetic Reference 
Field. 

11. The apparatus of claim 1 wherein said measurements 
of magnetic field data for three axes are measured for a 
preset time period. 

12. The apparatus of claim 5 wherein said measurements 
of magnetic field data for three axes are measured for a 
preset time period. 

13. A computer-implemented method for autonomous 
position determination of a vehicle comprising the steps of: 

(a) a program controlled processor inputting a measure- 
ment of three axes magnetic field data, said magnetic 
field data being measured for a preset time period; 

(b) said processor determining whether a sufficient sepa- 
ration exists between the magnetic field vector of the 
measurement and a direction of a geomagnetic dipole 
according to a dipole characterization of the geomag- 
netic field, to insure convergence of a correction 
scheme to a position solution, wherein the dipole 
characterization is a first order of a series of nth order 
spherical harmonics representing the geomagnetic field 
as a gradient of a scalar potential function; 

(c) upon a determination of non-existence of the sufficient 
separation, step (a) and step (b) being repeatedly per- 
formed until a sufficient separation exists between the 
magnetic field vector and direction of the dipole; 

(d) a program controlled processor calculating an inverted 
dipole solution comprising two position solutions, each 
position solution comprising a position vector for the 
vehicle and a vehicle distance from the center of the 
earth, for each measurement of three axes magnetic 
field data; 

(e) storing each position solution of the inverted dipole 
solution in a separate set in memory; 
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(f) repeating steps (a) through (e) for each measurement 

of the magnetic field data measured for the preset time 
period thereby creating two sets of position solutions 
propagating two position paths for the vehicle with 
different directions; and 5 

(g) said processor determining which of the position paths 
is the true position path of the vehicle. 

14. The method of claim 13, step (d) comprising the 
substeps: 

(dl) determining two unit position vectors of the two io 
position solutions in terms of the magnetic field 
vector and the direction of the dipole; 

(d2) said processor calculating three possible solutions 
for the vehicle distance from the center of the earth 
from the magnetic field vector, the direction of the 15 
dipole and a total dipole strength of the geomagnetic 
field; 

(d3) testing the three possible solutions to obtain one 
having a real nonnegative value as the vehicle dis- 
tance from the center of the earth; and 20 

(d4) said processor applying a correction scheme to 
each of the two position solutions of the inverted 
dipole solution until a stopping criteria is satisfied. 

15. The method of claim 14, step (g) comprising the 

substeps: 25 

(gl) determining a velocity vector for each position 
solution in each position path of the vehicle; 

(g2) for one position path, calculating a total energy for 
the position solution obtained from the first measure- 
ment of magnetic field data of the preset time period; 30 
(g3) for the same position path, calculating a total energy 
for the position solution obtained from the last mea- 
surement of magnetic field data of the preset time 
period; 35 

(g4) computing a difference of the total energies com- 
puted in steps (g2) and (g3) to obtain a first total energy 
difference; 

(g5) obtaining a second total energy difference by repeat- 
ing steps (g2) through (g4) for the other position path; 40 
(g6) testing whether one total energy difference is greater 
than a noise factor, a constant predetermined value; 

(g7) upon determining that one total energy difference is 
greater than the noise factor and one is less that the 
noise factor, selecting as the true position path of the 45 
vehicle, the position path having the total energy dif- 
ference less than the noise factor. 

16. The method of claim 14 wherein the correction 
scheme comprises a successive substitutions method having 

a stopping criteria of convergence within a number of 50 
iterations. 

17. The method of claim 16 wherein step (d4) comprises 
the substeps: 

(dl4) said processor inputting from memory Gaussian 
coefficients for all orders of the nth order model of 55 
series of spherical harmonics to iteratively calculate the 
higher order terms of the model with each iteration of 
position data. 

18. The method of claim 17 wherein said correction 
scheme further comprises following completion of said 60 
number of iterations for the successive substitutions method 
with applying a Newton-Raphson method for a second 
number of iterations. 

19. The method of claim 15 wherein the preset time period 

is no more than ten (10) minutes. 65 

20. The method of claim 19 wherein the preset time period 
is five minutes. 


21. The method of claim 17 wherein said number of 
iterations is 20. 

22. The method of claim 18 wherein said number of 
iterations for the successive substitutions is 20 and the 
second number of iterations is 5. 

23. An computer-implemented apparatus for autonomous 
position determination of a terrestial entity comprising: 

a magnetometer for providing a measurement of magnetic 
field data for three axes; 

an analog to digital converter for converting measure- 
ments of magnetic field data into computer readable 
form; 

a memory comprising a hemisphere location of the entity; 

said memory further comprising attitude data; 

a processor comprising a translation module translating 
the computer readable magnetic data from body coor- 
dinates to a fixed reference coordinate system by using 
the attitude data; 

said translated magnetic field data being stored in 
memory; 

said memory further comprising Gaussian coefficients of 
an nth order model of series of spherical harmonics 
representing the core portion of the magnetic field as a 
gradient of a scalar potential function; 

said processor comprising an inverted dipole solutions 
module for calculating an inverted dipole solution 
comprising two position solutions, each position solu- 
tion comprising a position vector for the entity and an 
entity distance from the center of the earth, for the 
measurement of three axes magnetic field data; and 

said processor comprising a position determination mod- 
ule for selecting a true position from the two position 
solutions based upon the hemisphere location. 

24. An computer-implemented method for autonomous 
position determination of a terrestial entity comprising: 

(a) obtaining a magnetic field measurement; 

(b) determining whether a sufficient separation exists 
between the magnetic field vector of the measurement 
and a direction of a geomagnetic dipole according to a 
dipole characterization of the geomagnetic field, to 
insure convergence of a correction scheme to an 
inverted dipole position solution, wherein the dipole 
characterization is a first order of a series of nth order 
spherical harmonics representing the geomagnetic field 
as a gradient of a scalar potential function; 

(c) upon a determination of existence of the sufficient 
separation, calculating an inverted dipole solution com- 
prising two position solutions approximately 180 
degrees apart, each position solution comprising a 
position vector for the entity and an entity distance 
from the center of the earth, for the measurement of 
three axes magnetic field data; and 

(d) selecting a true position from the two position solu- 
tions based upon a hemisphere location stored in a 
memory. 

25. An computer-implemented method for autonomous 
position determination of a terrestial entity comprising: 

(a) a magnetometer obtaining a magnetic field measure- 
ment in body coordinates; 

(b) said magnetic field measurement being converted to 
computer readable form by an analog to digital con- 
verter; 

(c) said magnetic field data being translated to a fixed 
reference coordinate system by using attitude data 
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stored in a memory, said memory further comprising a 
hemisphere location of the entity 

(d) determining whether a sufficient separation exists 
between the magnetic field vector of the measurement 
and a direction of a geomagnetic dipole according to a 
dipole characterization of the geomagnetic field, to 
insure convergence of a correction scheme to a position 
solution, wherein the dipole characterization is a first 
order of a series of nth order spherical harmonics 
representing the geomagnetic field as a gradient of a 10 
scalar potential function; 

(e) upon a determination of non-existence of the sufficient 
separation, step (a) through step (d) being repeatedly 
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performed until a sufficient separation exists between 
the magnetic field vector and direction of the dipole; 

(f) calculating an inverted dipole solution comprising two 
position solutions approximately 180 degrees apart, 
each position solution comprising a position vector for 
the entity and an entity distance from the center of the 
earth, for the measurement of three axes magnetic field 
data; and 

(g) selecting a true position from the two position solu- 
tions based upon the hemisphere location for the mea- 
surement. 



